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AERODYNAMIC INFLUENCE COEFFICIENT METHOD 


USING SINGULARITY SPLINES 

By J. E. Mercer, J. A. Weber 
and E. P. Lesferd 
The Boeing Company 


1.0 SUMMARY 


A new numerical lifting surface formulation, including 
computed results for planar wing cases is presented. This 
formulation, referred to as the "vortex spline" scheme, combines 
the adaptability to complex shapes offered by paneling schemes 
with the smoothness and accuracy of loading function methods. The 
formulation employs a continuous distribuution of singularity 
strength over a set of panels on a paneled wing. The basic 
distributions are independent, and each satisfies all the 
continuity conditions required of the final solution. These 
distributions are overlapped both spanwise and chordwise (termed 
"spline") . Boundary conditions are satisfied in a least square 
error sense over the surface using a finite summing technique to 
approximate the integral. The current formulation uses the 
elementary horseshoe vortex as the basic singularity and is 
therefore restricted to linearized potential flow. As part of the 
study, a non planar development was considered (i.e., interference 
shell body representation) , but the numerical evaluation of the 
lifting surface concept was restricted to planar configurations. 
Also, a second order sideslip analysis based on an asymptotic 
expansion was investigated using the singularity spline 
formulation. 



2.0 INTRODUCTION 


Various lifting surface schemes have been proposed for. the 
analysis of arbitrary wing planforms. Most of these methods treat 
this fluid dynamic problem by relating one flow parameter to 
another (e.g., downwash to pressure jump). A majority .of> the 
approaches can be characterized as constructing a solution from a 
composite of basic, assumed distributions of pressure, vorticity, 
or potential, all of which satisfy the governing differential 
equation. The boundary conditions of the problem determine the 
relative magnitudes of the various distributions., each of these 
assumed distributions can cover the entire lifting surface or 
merely portions or it, depending on the particular numerical 
scheme. 

Although many of the schemes arise from the above described 
common origin, in practice there can be great differences among 
them. As an example, consider the constant pressure panel (Ref. 
1) and pressure mode (Bef. 2) schemes. The constant pressure 
panel scheme uses a single type of distribution that is of 
constant pressure ranging over some small region of a wing. A 
solution using this method is one that consists of a set of steps 
in the pressure level which can only roughly approximate the true 
solution. However, the method is applicable to both subsonic and 
supersonic flows and ranges ever a wide variety of planform 
shapes. The pressure mode scheme, on the other hand, uses several 
different types of distributions ranging ever all or at least a 
large portion of a wing. The solution using this method can yield 
very accurate pressure results. However, the method requires 
careful attention to discontinuities in supersonic flow and it is 
restricted in the geometry cf planforms which it can treat. In 
choosing between these two schemes, trades between pressure 
details and generality of configurations must be weighed. 

The FLEXSTA8 system (Ref. 3) which performs the stability and 
control analysis of elastic arbitrary wing-body configurations was 
developed using the constant pressure panel scheme because of the 
schemes planform generality. Using this aerodynamic method, the 
program generates not only symmetric stability derivatives (e.g., 
C L„' C M„' etc.), but also antisymmetric stability derivatives 
(e?g., C| 0 , Cflg , etc.). It was soon recognized, however, that 
configurations with small wing dihedral could not be accurately 
analyzed for certain of their antisymmetric derivatives, 
specifically the rolling moment derivative, C| g . Upon 
examination, it was discovered that the accuracy of the constant 
pressure panel scheme was not adequate for the analysis of these 
configurations. In an attempt to establish the requirements 
necessary to treat such cases, a study was performed (Ref. 4). 
This study, based on an asymptotic expansion of the differential 
equation in terms of the governing flow variables (e.g., angle-of- 
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attack and angle of sideslip), revealed that gradients of the 
pressure distribution appear to play as active a role in the 
sideslip problem of configurations with small dihedral as the 
pressures themselves. This result confirmed the belief., that the 
constant pressure panel scheme could not treat this problem. 
Further investigation of other lifting surface schemes which might 
be applicable was initiated. This study concluded that no 
existing scheme could provide the needed pressure gradient 
accuracy required for the small dihedral configurations and at the 
same time, provide the desired generality. Accordingly, research 
to create a new schema was begun. The result of this study was 
the development of the "vortex spline" scheme: a method which 

satisfies the above mentioned requirements. The key features cf 
this method are: 

1) The use of independent, overlapping singularity distributions 
which are arranged such that the final solution for the 
pressure will automatically be continuous in both value and 
slope (unless specifically required not to be at certain 
locations) 

2) The use of wing paneling to provide a basis for the extent of 
coverage for each distribution, thereby giving it the 
generality of the constant pressure panel method 

3) The satisfaction of boundary conditions using a least square 
error technique 

4) The application to both subsonic and supersonic flow 

b) The use of "special ^ach lines" for paneling in supersonic 
flow. 

The current formulation of the scheme is restricted tc 
linearized, steady, potential flew. The restriction to linearized 
flow is due mainly to the use of the elementary horseshoe vortex 
as the basic singularity strength. This singularity was selected 
because of its automatic treatment of the linearized wake boundary 
condition and its compatability with the existing FLEXSTAB 
aerodynamic representation. , 

This report includes: 1) the derivation of the numerical 

scheme for flat wings of zero thickness along with computed 
results, 2) the description of a compatible interference shell 
representation (Hef. 1), 3) the description of a second order 

numerical scheme to compute interaction effects such as those 
which dominate for small dihedral, and 4) an alternate supersonic 
scheme studied initially (included in the appendix). 
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. 3 .0 NOMENCLATURE 


amplitude of quadratic -addition to panel edge, also 
used as a coefficient of the equation for a panel 
edge. -tr;r 

strength of Cauchy singularity 

area associated with ith control point 

aerodynamic influence coefficient, an elemerit of[A 

strength of logarithmic singularity 

strength of Mangier singularity 

strength of non-singular terms evaluated at y the 
singularity. 

weighted transpose of influence coefficient matrix. 

coefficient of the equation for a panel edge 

influence coefficient representing the effect of 
Subdivision I on Subdivision II $ :■ .. • 

coefficient of the equation for a panel edge 

average wing chord n . '• 

panel chord variation (x T£ - x TF ) - 

; m \ % £./' ■ : • 

section lift coefficient 
wing lift coefficient 
pressure coefficient 
square error 

integral of the kernel function in the x direction 
(linear in x variation ?of vorticity) 

•- -u ; 

integral of the kernel function in the x direction 
(constant in x variation in vorticity) 

numerical representation of Gj_ and G 2 


kernel function 



M 

N 

»o 

P 

r 

R J 


8 


U 

v 


W U 

W' 


w 


BC 


W r 


x 

X* 


y 

Y n 

t n 


Mach number, also, number of control points 
number of spline functions 
neutral point 
pressure 

cylindrical coordinate and also Mach line coordinate 

roots of Mach line - panel edge equation (intersections 
of Mach line with panel edge) 

area of integration and also Mach line coordinate 

strength associated with function , an element 
of{ 8 ) 

perturbation velocity in x direction divided by Ut© 
free stream velocity 

perturbation velocity in y direction divided by U,© 
perturbation velocity in z direction divided by U,© 
downwash on control point, i 

downwash on control point i, due to spline function, j 

downwash on downstream subdivision due to upstream 
subdivision (used in supersonic flow) 

downwash specified by boundary conditions 

downwash from spline functions 

Cartesian coordinate 

nondimensional x coordinate 

segment of chordwise vorticity spline variation 
(I is an index, not a power) 

Cartesian coordinate 

segment of spanwise vorticity spline variation 
(n is an index, not a power) 

(n is an index, not a power) 
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nondimensional y coordinate 
z Cartesian coordinate 

a angle of attack ■ 3; ' 

5 (subsonic) or ^M 1 -1 (supersonic) 

®1 ,C 2* C 3 sma H distance from singularity- '(used to evaluate G) 
T bound vorticity " 

Yj spline function vorticity distribution 

4P perturbation in pressure 

ACp jump in pressure coefficient across wing surface 

Aslope jump in panel edge slope ; V- 

Ax, % AX) quadratic additions to panel edges 
C running coordinate in x direction > 

n running coordinate in y direction 

o shed vorticity 

6 cylindrical coordinate 

• nondimensional 9 coordinate 

0 circumferential vorticity spline variation 

P density Of fluid 

♦ perturbation velocity potential 

♦d potential due to a point doublet 

• g potential due to a point source 

♦ v potential due to an elementary horseshoe vortex 

A$ jump in potential 
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Matrix Notation 


M 

K 1 

matrix of elements, Ajj 

.V * • • 

N 

diagonal matrix of elements, aj 


I s ) 

vector of unknown strengths, Sj 


w 

... vector of known boundary conditions 


[ f 

{ ) 

transpose of matrix 


vector 


[ r 

inverse of matrix 



Subscripts 


H 

homogeneous solution to differential 

equation 

i 

control point evaluation 


1 

lower surface 


L 

left panel edge 


LE 

leading edge 


ML 

Mach line 


0 

define location of singularity on shell surface 

P 

particular solution to differential 

equation 

R 

regular (non singular variable) and 
panel edge 

also right 

S 

singular variable 


TE 

trailing edge 


u 

upper surface 




4.0 REQUIREMENTS FOR NUMERICAL SCHEME 

Before developing a new scheme, a set of requirements and 
desirable features was established. The requirements were based 
on 1) ultimately incorporating the scheme into FLEXSTAB and' 2) the 
analysis performed in reference 4 (the latter because' the main' 
objective was to model flows around configuration with small 
geometric dihedral). The desirable features were based ' on 
experience with various lifting surface schemes. ! ** 

The requirement set. forth in reference 4 was that any 
numerical shceme that is to be used for interaction flow analysis 
(e. g. , sideslip-angle of attack) must be capable of predicting- 
gradients of pressure to the same degree of accuracy as most 
techniques currently predict pressure levels. This implies that 
the basic equivalent pressure distribution used to construct the 
solution must vary linearly as a minimum. A linear variation 
could approximate pressure gradients in the same way the constant 
pressure panel scheme approximates pressure levels. 


A further requirement, established by the ultimate desire to 
incorporate the final result into FLEXSTAB, was that the new 
technique must be adaptable to arbitrary planforms. ‘From 
experience with both pressure mode techniques and paneling 
techniques it was felt that treatment of arbitrary configurations 
can best be achieved through a paneling technique. Such 
techniques divide a configuration intc many parts (called panels) . 
A distribution of flow singularities is associated with each panel 
and therefore, these singularities conform to very complex 
geometries by mere adjustments to the panel sizes, shapes, and 
orient at ions . 


Experience using various numerical schemes also showed that, 
difficulties arise when attempting to model a lifting surface' 
which is in the wake of another (e.y., a tail behinu a wing ) i The 
source of the problem is singularities in the downwash field 
inherent in. various numerical formulations. Because of these 
singularities, the downwash influence is only correct at discrete 
points in the flow field. Fcr paneling techniques (which 
constitute the major source of this problem) this means that the 
trailing surface's paneling is determined by the leading surface's 
paneling. or vice versa (which ever provides the more stringent 
requirement). When dealing with arbitrary configurations many 
such paneling problems can arise thus causing much time to be 
spent establishing the paneling arrangement and in many cases, 
causing more panels to be used than would be for methods that dc 
not produce downwash singularities. This was not considered 
acceptable, therefore the requirement was established to eliminate 
as many of these singularities as practical. 
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As a desirable feature, it was felt that a basic scheme 
applicable to both subsonic and supersonic flew would be ideal. 
Although not a requirement, such an arrangement could reduce the 
coding needed to implement the scheme and at the same time allow 
a single mathematical model tc be analyzed all the way from 
incompressible flow through the supersonic regime. 

Finally, the most desirable general possible approach would 
be a single basic scheme which could be used for both analysis and 
design. This again would greatly reduce coding and at the same 
time, permit the user to almost arbitrarily switch analysis and 
design functions while designing an aircraft. As an example, a 
designer could generate a shape for a given flight condition and 
then analyze this configuration "off design" all in one run 
greatly speeding up the design process. 

Summarizing, then, the general requirements are: 

1) Capability to. model pressure gradients 

2) Arbitrary planfprm modeling 

3) Minimum downwash singularities in flow field 
The general desirable features are: 

1) Unified approach for .both subsonic and supersonic flow 

2) ... A single basic scheme for both design and analysis. 



Before a 

ny development of a new scheme was begun, a 

stu 

dy cf 

the 

characteristics 

of both subsonic a 

nd supersonic 

flow 

was 

und 

er taken to 

esta 

blish 

some of the 

more detailed 

nuie 

rical 

con 

si der at ions 

. For 

subsonic linearized 

flow the properties 

can 

be 

summarized 

as: 

' 





1) 

A s ingu la 

rity in the 

pressure cr 

bound vorticity 

at 

the 


leading 

edge 

which 

varies as the 

reciprocal of the square 


root of the dis 

tance 

from the edge - 

this causes the 

flo 

w tc 


turn abru 

ptly a 

t the 

leading edge. 




2) 

T he K ut. ta 

condi 

tion i 

it the trailing 

edge (vanishi 

ng 

bound 


v or tici ty 

or 

pressure) - this i 

s required to 

ma k 

e the 


analysis 

mathem 

atical 

Lly unique and 

physically ari 

ses 

from 


viscous e 

f f ects 

• 






For super 

sonic. 

linearized flow then 

e are several fe 

at ur 

6£ tC 

be 

considered: 
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1 ) 


Hyperbolic nature - a solution can proceed from fore to aft 
on a configuration. 

2) United regions of influence - disturbances are limited to 
the interior of Mach cones and ttndto concentrate most 
strongly along Mach cone boundaries*' 

3) Differing edge conditions - depending on the sweep of an edge 

with respect to a Mach line, the character of the solution 
differs: supersonic leading edges exhibit regular behavior 

whereas subsonic ones show a singular behavior.. Subsonic 
trailing edges require a Kutta condition, whereas supersonic 
ones do not. ' •. "■% 

4) Discontinuous nature of flow across "special Mach lines'* 
(those emanating from planfcra edge breaks - Fig. 1) - along 
special Mach lines the surface pressure and vorticity 
distribution can exhibit discontinuous beha vior in levels and 
for gradients. Within regions bounded by these, special Mach 
lines the flow exhibits a smooth regular behavior. / j 

Various schemes use different means to represent the above 
characteristics and in many cases the characteristics are only 
approximated; yet, these approximations yield adequate results. 
Therefore, it is not an absolute modeling requirement that a new 
numerical formulation possess the ability to exactly represent all 
of the above characteristics; it is only necessary to carefully 
study the impact of each in the evaluation of the scheme. 

The next section contains the analysis used to develop a new 
numerical scheme, the vortex spline, along with the computed 
results. The appendix contains another method for supersonic flow 
only, the doublet characteristic box, which was investigated but 
not fully developed. 






5.0 VORTEX SPLINE SCHEME 


5.1 General Approach 


As outlined in the in ttoduct ion , most lifting surface schemes 
model flow problems by relating one flow parameter to another. 
For analysis problems the known parameter is the vehicle shape 
while the unknown is the pressure distribution. Since both the 
known and unknown parameters can be expressed in terms of velocity 
and velocity related parameters the forthcoming analysis will be 
derived using these variables. 

Assuming linearized flow, the perturbation pressure can be 
expressed as: 


6P * ~W U (1) 

where 6u is the perturbation velocity in the free stream 
direction with U^, being the free stream velocity. Dividing 
Equation (1) by the dynamic pressure, ^ » gives: 

6P ,_ 2 « H _ ( 2 ) 

I foo OO 00 

or 


c p = - 2 u 


(3) 


where Cp is the pressure coefficient and u = 6u /Va> is a 
nond imensional perturbation velocity. 


Across a lifting 
corresponding to a 
jump can be expressed 
parameter - the bound 


surface there is a jump in pressure, ACp t 
jump in the perturbation velocity, u. This 
in terms of another velocity related flow 
vorticity, y : 


y = 


- u- 


(4) 


where u u is the upper surface velocity perturbation and u. is 
the lower surface velocity perturbation. 
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Since Y corresponds to an antisymmetric shear 
lifting surface (Ref. 5, pg. 129) 


ac ross 


the 


therefore: 


i -- 



a rid : 


u 


u 



u 



( 5 ) 


( 6 ) 


AC p ■ 2(u u - Uj ) ■ 4u “ 2Y 


(7) 


From Equation {7) it is clear that the jump in pressure 
across a surface, assuming linearized ©flow , is directly related to 
the shear flow parameter - bound vorticity, Y . 

on a two-dimensional planar surface, there are two components 
of velocity which are discontinuous, u and v; there are 
corres pondingly two components of vorticity - Y corresponding to 
u and o corresponding to v {see Ref. 6, py. 221): 


o i -(v u - Vj ) - -2v 

whore' v is the nond imensiona 1 , perturbation velocity normal to the 
free stream. 


In treating potential flow, a and y are not independent 
and can be related through Helmholtz theorem (Ref. 6, pg. 221) as 
(Fig. 2) : 

5 


i? + 21 

BX ay 



(9) 


An elementary 
equa t ion : 


solution to the governing differential 


( 1 “ M J ) 





0 


( 10 ) 
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which produces shear flow corresponding 
obtained by integrating the basic doublet 
termed the "elementary horseshoe vortex" 
the linearized wake boundary condition of 
no jump in pressure. 


to vorticity can be 
(Ref. 5). This solution 
automatically satisfies 
a jump in potential but 


Using the elementary horseshoe cortex as the basic 
singularity, planar flow problems can be modeled as follows: 


w(x,y) - ff m. n ) K (C-x, n - y ) dC dn (11) 

8 

where w(x,y) is the known nondimens iona 1 downwash velocity (for 
small angles w is approximately a , the angle of attack) , x and 
y are field points, s is the region over which the singularities 
are distributed, T(C,n) is the unknown bound vorticity 

distribution, and K(c-x, n-y) is the horseshoe vortex kernel 
f unc ti on . 

For subsonic flow, the region of influence, s , and the 
entire singularity surface coincide. The kernel function for this 
flow regime is (Ref. 5, pp. 87 and 149): 


K(C - x, n-y ) - 


where 0 Z =1-N 2 . 


1 1 

C - x 

4* (n- yF 

^(C-x) 5 + * J (n - y) k 


( 12 ) 


For supersonic flow, the region of influence, s, is the 
portion of the singularity surface contained in the forward Mach 
cone emanating from the field point (x,y). The kernel function 
for this flow regime is derived in the same manner as was done for 
subsonic flow (Ref. 5, pg. 87) except that the analysis begins 
with the supersonic source (Ref. 5, pg. 88). Since the supersonic 
singularities can only influence locations within the aft Mach 
cone, the integration along the line of doublets does not range 
from x to oo but rather from x to the intersection of the line 

of doublets with the forward Mach cone -from the influenced point 

( i. o. , x - fi yfr'-y j? for a doublet line in the z=0 plane). The 

resulting kernel function is: 


0 


2 


K (c - x , n - y ) = 

is defined as K 


1 J f - x 

" 2 " 

*-*1 for supersonic flow. 


(13) 
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where 


As mentioned earlier, Y(c,n) # the unknown, c ( an- be 

represented in an approximate way as a series ot basic 
distr ibutions r j U , n ) , so that : 

' m . ’ “ 

Y (<,n) = ESj Yj (C , n ) "(14) 

J*l- . ■ 

where Sj is the strength of the distribution. This technique 

permits the integration tc be performed using the , known 

distributions of y.(c , n ) ; thereby reducing the problem ,tc an 
algebraic solution^ for Sj . The algebraic formulation is 
established by specifying the downwash at M>N "control points" 
(i.e., collocation points): 


"(x^yj) = wp i-1, ....M (I 5 ) 

where w A f s the downwash at x i ,y i . Substituting Equation (14) into 
Equation (11) and integrating yields: 

N 

w i = ^-,7 •••• M (16) 

a V ^ 

The aerodynamic influence coefficient, Apj = 3Sj , represents 
the integral in Equation (11) (with Sj =1) and gives J t.he downwash 
at point i due to unit strength of the jtbdistribution of 
vorticity. Equations (16) are a set of algebraic equation fcr Sj . 
For K equal to N, Equations (1b) are a determinate set of 
algebraic expressions with’ a unique solution for the S j . For K 
greater than N there are more equations than unknowns and, as will 
be shown later, the Sj can be determined in a least square error 
sense. Once tne Sj are computed the pressures and pressure jumps 
can be computed from Equation (7). 

Lifting surface schemes using the horseshoe vortex kernel, 
differ in the form of the basic distribution, Yj (t,n) that they 
assume. They also differ in the manner in which the boundary 
conditions are satisfied, through the location of the field points 
and/or the number used in relation to the number of unknowns. The 
next section describes the type of vorticity assumed for the 
vortex spline scheme. 
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5.2 Selection of Vorticity Distribution 


o Two-Dimensional Study 

Defore attempting a full three-dimensional study, 
various numerical techniques were first developed in two- 
dimensions. These two-dimensional studies supplied insight 
while at the same’ time providing a fundamentally simpler 
means to obtain numerical results than would be required for 
the full three-dimensional problem. The initial 
■ investigation concerned various vorticity distributions. In 
two dimensions, the basic vbrticity kernel function for 
downw ash is: 


K (x) 


1 1 
2 ir x -C 


therefore the downwash is: 



(17) 


(18) 


A study was ma.de using this expression as a two- 
dimensional formulation of the constant pressure panel 
scheme. With Y ( C ) = 1 for x^ < C < x 2 the integration can 
then be performed yielding: 



w(x) 


1 

2w 


In 


x - 

X - Xg 


(19) 


Figure 3 shows how the scheme is used to model the 
loading on a flat plate at angle of attack. 

The vorticity distribution is a series of step 
functions. The influence of the j^h-ste'p function, having 
strength, Sj , on the i th control point can be expressed from 
Equation (19) as: 







where 

Aij = -i- i n |fi, X J I 
2 " *i - x J+1 l 


( 22 ) 


The location of the control points is computationally 
determined so that the numerical results will be correct over 
a wide range of problems. Since there is one control point 
and one value of Sj corresponding to each panel. Equations 
(16) constitute a determinate set of linear algebraic 
equations (tt=N) which can easily be solved for the strengths 

S J * 


The thing to note about the constant pressure panel 
scheme is that the downwash field for the panels contain a 
logarithmic singularity at each panel edge. This means that 
the control point placement on the panel is critical since 
the gradients of the downwash can be large. 

Figure 4 shows the results of the flat plate analysis 
for subsonic flow. These results show the large downwash 
gradients and also show that there is no way to specifically 
invoke the Kutta condition at the trailing edge. This 
condition is implied through the control point location. 

In spite of these problems the results are fairly good; 
but, however, they are not sufficient for the interaction 
problem (as pointed out in the introduction). 

In order to improve upon these results it was apparent 
that two things had to be done - 1) remove as many of these 
logarithmic singularities as practical, and 2) establish some 
definite means to introduce the Kutta condition in subsonic 
flow. The first item was achieved by noting that the 
logarithmic singularities arise from jumps in the vorticity 
strength therefore a new distribution had to be formulated 
which had continuity of vorticity level. This implied that 
the vorticity had to posess at least a linear variation 
(consistent with the requirements set forth in Ref. 4, and 
outlined in the Introduction). There are two means to obtain 
a continuous distribution using a linear distribution - one 
is to place separate continuity constraints on the solution 
across panel edges using a method such as Lagrangian 
multipliers; the other is to form the functions so that they 
automatically provide continuous distributions. This second 
means was selected because of the inherent simplicity in the 
mathematical formulation. Figure 5 shows the basic function 
used for the new two-dimensional numerical scheme. . This 
figure also implies that there is still a jump in the 


19 



FIGURE 4 - TWO DIMENSIONAL CONSTANT PRESSURE PANEL RESULTS 
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FIGURE 5 - LINEAR VORTICITY REPRESENTATION USED CHORDWISE 








vorticity at the leading edge which causes a logarithmic 
singularity. The only way to remove this singularity at the 
leading edge is to use the proper asymptotic form which is 
related to the reciprical of the square root of the distance 
from the edge (Figure 4). The new scheme was tested with and 
without the proper singular form, the conclusion was that 
there was sufficient accuracy for lift and moment without 
this singular variation in the vorticity. 

Figure 5 shows that the triangular distribution of 
vorticity provides a natural means of introducing the Kutta 
condition by omitting the last function. This emission 
forces the vorticity to vanish at the trailing edge. 

The vortex spline formulation is based on a different 
means of satisfying the boundary conditions which will be 
discussed later; but the fundamental approach is to compute 
the downwash at discrete locations using the basic triangular 
distributions of vorticity, Y , in Equation (18) - with the 
limits changed to include two panels per distribution (Figure 
5), Figure 6 shows the results for subsonic flow. 

O Basic Three-Dimensional Function 


For the three dimensional formulation 
components of vorticity must be considered, 
pointed out earlier, these two components a 
in potential flow. This means that it 
arbitraily legislate that both components 
all directions. Using Equation (9) to solv 
of y gives; 
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Summarizing then: 
y - linear in x, quadratic in y 
o - linear in y, quadratic in x 

(A plot of o versus y would be analogous to Figure 5- 
a series of linear functions.) 

Hith this above variation the only places where 
singularities would occur on a simple rectangular wing would 
be at the leading edge (due to a jump in T ) and at the tips 
(due to a jump in a ) for subsonic flow and at the tips only 
for supersonic flow (Figure 7) . 

Since the elementary horseshoe vortex is best described 
in terms of y , the criteria established above for this 
component of vorticity was the one used as the basis to 
derive the three dimensional distribution function presented 
next. 

For the x variation of 7 in three-dimensions, the 
requirement shown above is the same as that for two- 
dimensions. However, the y variation has to be quadratic or, 
in more fundamental terms, continuous in both value and slope 
in order for a to be continuous in value. The logical 
means to form such a surface distribution seemed to be as a 
product of two functions - one providing the desired linear 
x variation and another providing the desired quadratic y 
variation. To see how this y variation was formed, consider 
the basic triangular form used in two-dimensions. This 
distribution has, in general, three discontinuities in slope. 
By adding three discontinuous quadratic distributions to this 
basic form, these discontinuities can be removed. Figure 8 
shows how this is accomplished. The requirement for each of 
the quadratic functions is that they begin and end at the 
centers of panels (spanwise) with zero value and slope, and 
also have a discontinuity in slope at the panel edge which 
cancels the discontinuity in slope of the triangular 
function. The overall function then becomes six separate 
quadratics which cover four panels spanwise. The equations 
of these quadratics can be computed from the imposed boundary 
conditions cn each segment. Table 1 gives these results. 
(As a note, these spanwise functions could have also been 
developed using three quadratics over three panels specifying 
the slope and value matches at the two inner panel edges, and 
also requiring the distribution to vanish in both value and 
slope at the outside edges. The additional requirement of 
unit maximum value would then completely specify the 
distrib ution. ) 
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Note: no downwash singularity at leading edge, if edge is supersonic 


FIGURE 7 - LOCATION OF DOWNWASH SINGULARITIES WITH VORTEX SPLINE 
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FIGURE 8 - QUADRATIC DISTRIBUTION FORMED FROM LINEAR FUNCTION AND QUADRATIC ADDITIONS 

































If X^is used to represent the chordwise variation of Yjj 
and Yj is used to represent the spanwise variation of , the 
distribution of vorticity could be expressed as a aoubly 
subscripted variable, Y^, j . However, in order to conform 
with the singly subscripted Y^ of Equation (14), a unique 
relationship between k and i,j was assigned so that: 


\ = Y i’J " X i Y J 


. . N 


where Yj consists of the six quadratic segments Yj n (n = 
1,2,.. .,6) shown in Table 1 and Xj. consists of the two linear 
segments X® (m = 1,2) shewn in Figure 5. Thus, the analytic 
description of Y^ = Yj^j consists of twelve discrete 
parts, X® Xj n (m = 1,2; n = 1,2,..., 6). These parts are 
distributed over eight panels (four spanwise, two chordwise) 
as shown in Figure 9. 


For a set of rectangular panels, is only a function of 
the streamwise coordinate, x. By definition, Yj is only a 
function of y, providing the spanwise edges are aligned with 
the free stream (which by the panel definitions, will always 
be the case). If the panels are not rectangular, Xj becomes 
a function of both x and y so that it varies as the fraction 
of local chord: 


In the above equations, x. 


leading 


trailing edge coordinates of the panel covered by segment 
(X^ ) ; this is illustrated in Figure 5. 


Except 
functions of 


rectangular panels, X ^ and 
If the fraction of local chord 


x - X LE 
*TE ‘ x LE 


then Xj. becomes 
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FIGURE 9 - SURFACE SPLINE DISTRIBUTION OF VORTICITY 
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(30) 


The above representation holds for all leading and 
trailing, panel edge shapes. 

The bound vorticity distribution defined by Equation 
(25) (termed spline) and shown in Figure 9, produces no 
singularities in the dovnwash. This function is scaled up 
and down in the solution exactly the same way as the constant 
pressure panel scheme usinq a multiplicative constant which 
is the unknown. 

The next sections describe how a wing is subdivided into 
panels and how the basic splines are overlapped using these 
subdivisions. 

o Modeling Vorticity on a wing Using Vorticity Splines 

Figure 10 shows a paneled wing with both a spanwise cut 
and a chordwise cut to display the chordwise and spanwise 
overlapping of the basic spline functions (assuming unit 
amplitude). For the five chordwise functions and the six 
spanwise functions shown, there will be a total of thirty 
(N=30) spline functions distributed over twenty-five wing 
panels. The components of a single function are identified 
(shaded area) to depict its relationship to the overall 
arrangement. This figure also shows that the Kutta condition 
has been imposed (zero value of function at trailing edge). 

As shown in Table 2, special functions are used at the 
root and tip, the functions for symmetric flow will be 
discussed first. 

The first root function (root function 1) is the right 
half (spanwise) of a normal function. The second root 
function (root function 2) is composed of segments 
y3 y** Y 5 and yk of a normal function with y^ of the 
corresponding function on the left half of the wing added to 
y2 of the right side function. Adding these two terms 
accounts for the overlapping of the right and left halves. 
(This can be done since the magnitude of functicr. 
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Root Function -• Symmetric 
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(y ranges from 0 to 1 
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Note: The normal function represented by segments 4 through 6 at the root 


meets the proper symmetry requirement of vanishing slope. 

Root Function 1 - Antisymmetric 

Segments 1 through 6 do not exist 
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TABLE 2 - SPECIAL SPANWISE FUNCTIONS 


































Root Function 2 - Antisymmetric 
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TABLE 2 - CONTINUED 
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Tip Function - 2 

























corresponding to the segment is the same as that 

corresponding to the segment-by symmetry. The reason for 

this treatment will be explained in the next section.) 

The next to last tip function (tip function 1) has Y^ - , 
Y^ , and y3 as normal segments. y6 does not exist; y^ and 
y5 are one continuous guadratic function, the boundary 
conditions of which are: 1) continuous slope and value at 

left edge of the tip panel, and 2) zero value at the tip. 

The last tip function (tip function 2) has Y^ as a 
normal segment, but y2 and y3 , the only remaining segments, 
are determined such that when this tip function is added to 
the other functions, the combination will produce a constant 
value of y it all the strengths (multiplicative 

coefficients) are the same. This requirement was set 
because, in supersonic flow, the vorticity strength can be 
constant approaching a pointed tip. This tip distribution 
modeling capability also tended to yield the best numerical 
results of any of the variations studied (for both subsonic 
and supersonic flows). 

For antisymmetric flow (e.g., rolling) the root 
functions have to be altered. Figure 11 shows these modified 
spanwise variations. The first root function (Y^) has been 
deleted. The most inboard function (Y 2 ) , for the 

antisymmetric case, is the same as a reflection of the next 
to last tip function (Y^). The remaining functions across 
the span are the same as for the symmetric case. 

Table 2 gives all the coefficients for the special 
functions just described. 

Figure 12 shows an example of how two functions are 
overlapped - that is share common panels. This illustration 
depicts one of the many possible ways that overlapping can 
occur. 

Arbitrary flow problems involving symmetric 

configurations can be decomposed into a symmetric and an 
antisymmetric flow problem. In this way, only one half of 
the wing need be treated for each of these reduced problems. 
For symmetric flow, the unknown strengths have the same 
values for corresponding functions cn opposite sides of the 
configuration's plane of symmetry. For antisymmetric flow, 
the strengths have the same magnitude but opposite sign for 
corresponding functions. Using these properties to formulate 
the problem results in the concept of real and image 
relations and the satisfaction of boundary conditions on one 
side of the plane of symmetry only (the other side 
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FIGURE 11 - ANTISYMMETRIC SPANWISE SPLINE DISTRIBUTION 
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FIGURE 12 - OVERLAPPING Of SURFACE DISTRIBUTIONS 



automatically satisfied by the appropriate symmetry 
relation) . The symmetric problem is constructed as the sum 
of the influences from the real and image functions, where 
real and image implies corresponding functions on opposite 
sides of the symmetry plane. Tho antisymmetric problem is 
constructed as the difference of the influences from the real 
and image functions. This scheme reduces the number of 
calculations to be performed by half. 

o Wing Paneling 


As mentioned earlier, the basic paneling procedure is to 
cut the wing streamvise to form panels so that two panel 
edges are parallel to the free stream. This arrangement 
allows the kernel function integration to be performed more 
readily by introducing simple limits; it also simplifies the 
expression for the Y variation of . the basic vorticity 
distribution (making it a function of y only). 


The initial assumption for paneling was that the new 
technique use the same trapezoidal paneling as the constant 
pressure panel method (Ref. 1) . Quickly, it was realized 
that any non-rectangular planforas would have discontinuities 
in the slopes of the fore and aft panel edges. Such 
discontinuities would introduce logarithmic singularities in 
the downwash field, including the wake. The source of these 
singularities can be seen by inserting the expression for the 
Y spline distribution (Eg. 25 ) into Equation ( 23 ) • 

xj^y),* 

o(x, y) m ~h f xU, y)*(y) d c (31) 

XLE(y) 

The upper limit of the above integral depends on whether 
o is computed on the distribution (limit equal x) or behind 
the distribution (limit equal Xfg(y) ) . 


The above equation shows that if any of the variables 
T(y)» X(x,y)» Xle (y) or XTE(y) has a discontinuity in slope with 
respect to y then o will be discontinuous in value 
producing the logarithmic singularity. Of course Y(y) is 
generated so that it is continuous in both value and slope so 
will not be of concern. This shows that planform * breaks, 
which cause jumps in panel, slopes, will produce the 
undesirable downwash singularities. 


A method 
eliminating 
procedure and 
edges. The 


was developed to eliminate this problem by 
the planform breaks. Figure 13 shows the 
gives the resulting equations for the planform 
procedure shown is exactly the same one used to 
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FIGURE 13 - QUADRATIC PANEL REPRESENTATION ALONG A PLANFORM BREAK 



change the basic, linear distribution of vorticity into the 
quadratic distribution (Figure 8). 

An additional benefit of going to the quadratic edge 
panels is that they provide a more general modeling element 
to handle arbitrary planforms. Later results in this report 
will show their usefulness in modeling a circular wing. 

The basic panel then is one with two parallel, 

streamwise edges with two quadratic edges; the quadratic 
changes coefficients at mid span. With this paneling 

arrangement there will be no downwash singularities in the 
wake (except from the tip) ) and therefore, meets the desired 
requirement. (The quadratic paneling could also have been 
developed using one continuous quadratic across each panel 
edge in order to be compatible with the three spanwise panel 
developments mentioned earlier.) 

The paneling described above is the one always used in 
subsonic flow. For supersonic flow, however, there are two 
choices of paneling - the one just described and another 
termed "special Mach line paneling." In this latter 

procedure certain Mach lines are identified, across which the 
vorticity level and/or gradient can be discontinuous (see 
Figure 1) . These lines eminate from planform breaks, tips, 
etc. Along these lines continuity conditions can be relaxed 
and downstream sections can be treated as separate wings in 
the wake of the upstream sections. Interior to each section, 
the paneling is treated as though that section were a 
separate subsonic wing (normal quadratic panels) , so that the 
Tj distributions for each subdivision do not overlap panels 
of an adjacent subdivision. 

Figure 14 shows a square wing at M equal V~2 comparing 
the normal geometric paneling with the special Mach line 
paneling. Note that the planform breaks are eliminated even 
when they are due to the special Mach lines. Also note that 
since there are no singularities in the wakes emanating from 
each subdivision, the paneling of the subdivisions can be 
independent (except at planform breaks where the edges have 
to be aligned to ensure that the quadratic edges fore and aft 
coincide) . 

Numerical results computed with the special Mach line 
paneling (presented later in this report) show that this 
procedure greatly improves the details of the pressure 
distribution and at the same time allows for even greater 
variation in paneling options then afforded by the normal 
geometric paneling. 
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5.3 Evaluation of Kernel Punction Integrals 

Equation (11) gave the integral expression which relates the 1 
downwash, w, to the vorticity distribution, y • Equations (12) 
and (13) gave the kernel functions appearing in the integral for 
subsonic and supersonic flow, respectively. - This section 
discusses the evaluation of the integral for both flow regimes and 
the additional problem of establishing the limits for the 
supersonic integral. 

Substituting the expression for Yj , the basic spline 
distribution of vorticity (Equation (25) ) , into the downwash 
integral . (Equation (11) ) gives the downwash at point i due to unit 
amplitude of the jth spline distribution: 


y R x TE(n) 

'ij * / / 

y L x LE(n) 


[x ( C , n)Y(n)j.K(C- * t , n - y^dCdn 



where [ XY ] . denotes the right hand side of Equation (25). Y L and 
*R are the reft and- right edges, respectively, of the panel 

associated with X«* Yj m (components of [XY]). Xje ( D ) and ( n ) 

are the trailing ana leading edges, respectively, of that same 

panel. The symbol f refers to the principle value of the integral 
(Hef. 5, pg. 132)r. The total downwash at x A , y^ is then the 
sum of the influences from all the integrals for the various 
spline distributions, each multiplied by its respective strength, 
Sj . in Equation (32), corresponds to Ajj in Eguation (16): 

»i • £ a ij 8 j (16) 

J"1 

Examination of the integral in Equation (32) shows that since 
yr and y^ are constants, due to the fact that they are parallel 
edges, the integration can best be performed on C first. 

Returning to Equations (26) and (27) for X, it is noted that both 
can be written as: 


Y . + c-f ( n ) 
c ( n ) 


<3 3) 
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where the variables of integration, 
( y, respectively, 
replacing Ajj for 


U 


and 


2 


C and n t replace x and 
Substituting Eguation (33) into (32) and 
yields two general forms: 

y R x TE(n) 

± f CK(t-x 1 »n-y l )dC (34) 

y L x LE(n) 

5* tin) ? E < n > 

/ y( ">c7kW K( '” 


n - y t )dc 


[ LE(n) 


(35) 


For arbitrary edges the integrations indicated by Equation 
(34) and (35) cannot be performed analytically for both { an 
o . The integration in the C direction can be performed 
analytically, however, leaving the integration in the n 
direction to be performed numerically. Even for the special cases 
when the integration can be performed entirely analytically (e.g., 
rectangular panels) , it is faster to perform the cuter, r> 
integration numerically. This is because of the number of 
hyperbolic functions, etc. which arise from the analytical 

expression and must be evaluated by the computer. 


Letting : 



G,(n) 



J 

X LE 


t KU- *i » n - Vi) 

(n) 


dC 


and 


x TE^ n ^ 

C 2 (n) - / K(C- , n-y^dC 
x LE( n ) 


(36) 


(37) 


the resulting expressions for G(n) in subsonic flow can be written 
as: 1 C [ r2 (x.+ c) , 3 

G l( n ) " ~u ( y* - n) 2 [ V ' 2 \J(x 4 - c) 2 + fl 2 (y 1 -n) 2 J 


G 2( n >- 4T~ (y^-T) 2 ^ * *1 ’ c ) 2 8 2 <V n > 


- 1 n|( x |-C)^(x 1 -c) 2 + fi 2 (y r n) 2 l c [ x 

t 3 x TE^ n ) 

O X LE ( n ) 


t* X TE^ ^ ^ 


le(°) 


(38) 


(39) 
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and in supersonic flow as: 


■ -rr 


xi + ; 


e f yi • n p 


. yin 


f 2 2 2 

X\ - O*-, (*1-0 -8 (y 1 -n) 


*• m x TE^ n ) or x ML^ n ^ 

C ’ x le(^>) (uo) 

C= x te f ^ ) or x ML^ n ^ 


2(n) = " ^ x l-o 2 - fi 2 (y t -n) 2 


S * x LE^ r| ) 


(41) 


The notation at the right side of the above expressions is 
the standard one for integrals, indicating that the expression is 
first evaluated with c replaced by the upper limit and then, 
subtracted from that, the expression with t replaced by the 
lower limit. 

Under certain circumstances Equations (38) through (41) can 
become singular due to the (y^ -n) 2 term in the denominator. The 
order and presence of a singularity in the integrand depend on the 
relationship of y, to the two side edges of a panel, y L and y R , 
and the relationship of Xj to the other two edges *T.refo) and xte(o)» It 
is obvious that if the control point y value falls to the left of 
y L (the left edge) or to the right of y R (the right edge) the 
integrals involving (y.-n) will never be singular. For this case, 
the numerical integration with respect to n from y L to y R can 
be performed routinely using Gaussian quadrature. It is also 
true, but less obvious, that if the x value of the control point, 
Xj , lies upstream of (yj ) then there will be no singularity 
regardless of the y t value (this can be verified by taking the 
limit of the expression evaluated at both limits as n approaches 
y, ). For this last case, however, care must be taken that none 
or the quadrature points falls too close to y t , when y A falls 
between y L and y R . Such proximity could give rise to numerical 
problems. Other than taking this precautionary measure, the 
integration is routine. 


If the y. value falls between y L and y R and the x* value is 
downstream of the leading edge,x L g(y| ), then the G ( n ) functions 
cannot be integrated directly. ror this case there will be three 
basic singular forms: 1) a logarithmic singularity (for only the 
linearly varying terms, Gi), 2) a Cauchy singularity (1/yj-n fcrm) , 
and 3) a Mangier singularity (l/(yi -n ) 2 form). The first 
singularity is integrable in a Riemanian sense,, but the latter two 
are not. These singular forms arise from taking the limiting form 
of the kernel function as the point of evaluation approaches the 
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(plane of the singularity. There would be no singularities if the 
■ limiting process had been taken after the integration. The 
'advantage of taking this limit before integrating is the resulting 
simpler expression. The cost of taking this limit is the 
introduction of the singular forms. Although the Rangier and 
Cauchy expressions cannot be numerically integrated, they can be 
eliminated from the integrand and treated separately (Ref. 7). 
There are several ways to accomplish this, but basically all the 
methods assume that the integrand can be expanded about the 
singularity as: 


‘M 




+ A L ln 


(y t -n)U A 0 4 A 1 (y l -n)+- Ag^-n) 2 * 


( 42 ) 


Here (T ( n ) represents either G x (n) or Gg ( 0 ) multiplied by 
the remaining if functionality (from Equations (34) and (35)). 

Expansion 42 can be performed either analytically or 
numerically. The analytical method is covered in Reference 7, 
so it will not be repeated here. 

Numerically “( n ) can be expressed as: 

G(n) « G s (n) ♦ G R (n) (43) 


where Gg ( n ) contains all the singular forms and Gp ( n ) is 
regular. In the vicinity of the singularity, the leading terms of 
the series in Equation (42) dominate, therefore: 


+ A 0 (44) 


The only term of the above series which is antisymmetric 
about the singularity is the Cauchy term, Ap/yj-n , therefore, 
it is possible to evaluate the coefficient of this term 
independently by computing G(n ) at two points clcse to the 
singularity" n»y^ + e^and n = yy€^: 

A c 2 - G(yi + e x) ] (45) 

The choice of is a function of the number of significant 

figures the machine uses to store numbers and, to a lesser extent, 
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A l la 


A 

yi'-n 

(In region of singularity) 


V 



the manner in which the expression for G(n ) is constructed. (The 
accuracy of G( n ) near the singularity is best if the origin of 
the axis system is located at the singularity. This origin 
location avoids the numerical problem of adding a small number to 
a large one.) 

Since the remaining terms of Equation (44) are all symmetric 
about the singularity, they must be evaluated simultaneously. The 
three unknowns Am ' A L # an( * A 0 can be numerically computed by 
evaluating ”(n ) at three points. Because the Cauchy term has a 
rather large influence, numerical error can be reduced by 
eliminating this antisymmetric term (as well as the__ other 
antisymmetric terms) from the calculation by evaluating G(n ) at 
equal distances on both sides cf the singularity and summing the 
influences. This gives three equations for the three unknowns: 


2 

f 2 A M 
1 


+ 2 In I 


£i| a L 


+ 2A 0 - Cty* + Ci) + G (*! " e l ) (46a) 


T a m 

e 2 


2ln|c 2 |A L + 2A 0 = G(y 1 + e 2 ) + G(yi-e 2 ) 


(46b) 


2. 

"p A M 

C 3 


+ 2 lnl Eol A 


3 1 a L 


2Aq 


g (y t + e 3 ) + 


e(» t - e s 


(46c) 


These equations can be solved for and A^ • Aq will not be 
required as will be shown next. 

Using the above results, Gg (n)r Equation (43), can be 
expressed as: 


G S (n) 


V 

(y 1 * r ») s 





- n 


(47) 


The above expression can be integrated directly using the 
integrals given for the Mangier and Cauchy forms in Reference 7. 


/ 

*L ( y i~ n ) 


2 dn = A m 


(TT-TrT ' (yr»L) 


(» 8 ) 




^ R (n) - G(n)- G s (n) - G <*) - 2 " 7 ^“ ‘ A L *» I V n I (51) 

The above expression picks up the Aq and remaining terms __which 
were dropped to compute the Aj^ , A q , and A^ . Gr (n) 
contains no forms which cannot be integrated numerically^ using 
Gaussian quadrature. Therefore, numerical integration of Gr ( rj ) § 
when combined with the results shewn in Equations (48) through 
(50), constitute the entire integration of G ( n ) to evalua te Ajj. 
(It is important to mention again that the allowable magnitudes 
for , Cg , and Cj depend on the number of significant figures 
carried for G(n) and are therefore computing machine dependent.) 

The above procedure was compared to a method which used the 
non-singular, out of plane expression for the kernel. This 
alternate procedure integrates the out of plane kernel using a 
control point very close to the singularity plane. For the same 
accuracy, the planar computational scheme, just described, was. 
about ten times faster. Both schemes used a self checking 
integrator which breaks the interval of integration into parts 
small enough to be handled by a four point Gaussian quadrature 
using eight points as a check (Reference 8). 

The C integration in supersonic flow is not always between 
the leading and trailing edges. Figure 15 shows the forward Mach 
lines emanating from the point * 1 , . These lines are shown 
cutting a panel, the significance of this cut is that portions of 
the panel ahead of the lines influence the point x A ,y^ , while 
those aft do not. In terms of the kernel function, the regions 
aft of the lines cause the values in the radicals to become 
negative making the expression for the integrand complex'. 
Physically, this means that the influencing region is limited to 
the area within the forward Mach lines, therefore, the limits of 
integration should be adjusted to include only this region. 
Figure 15 shows two distinct regions of integration, one from the 
panel leading edges to the panel trailing edges - Region I; and 
one from the panel leading edges to the Mach lines - Region II. 
The integration over a panel is the sum of the integrals over the 
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two types of regions (as will be shown in Figure 16, there can be 
several of each type on a single panel). 

The integration of the supersonic kernel in the C direction 
from leading edge to trailing edge was given in Equations (40) and 
(41). The integral from leading edge to Mach line is the same 
except for the upper limit, Xq^(n)» replaced by the equation of 
the forward Mach line: “ 


C ML " x i + *(n - /j ) 


— right going Mach line 

(52) 

+ left going Mach line 


This upper limit on C simplifies the expressions in Equations 
(40) and (4 1) by- forcing the radicals to vanish. The integrands 
for the n integration which result from this changing of limits 
one are therefore no more complicated to integrate than the ones 
already discussed^ The problem then is identifying the various 
regions on the panel so that the limits for the n integrations 
can be determined. 

For curved panel edges, a Mach line 
ways. Figure 16 shows some of these 
intersection characteristics can be noted 

1. After an intersection with a trailing ."edge’* the Mach line can 
intersect either the trailing edge again or the leading edge. 

2. After an intersection with a leading edge the Mach line can 
intersect either the leading edge again or the trailing edge. 


can intersect in various 
combinations. Several 
from the figure : 


3. Proceeding left to right (positive n or y direction) each 
intersection changes the region - from a Region I to a Region 
II from a Region II to a Region I and from a Region II to a 
Region 0 (Region 0 does not influence the control point) ; but 
there can never be a jump from a Region I to a -Region 0 or 
vice versa. * 


The above characteristics can be used to form a logical 
scheme to. identify the various regions and find the limits of 
integration. ^Usirig the quadratic equation for the planfcrm edges 
from Figure 13 and the equation for the Mach lines (52) expressed 
in a common coordinate system, a quadratic expression for the 
intersections can be obtained: 

,2 s o (53) 


x i + *(n - *i) - a - bn 
where a +- bn + cn^ is the 


cn’ 


of the panel edge 


equation 






The four roots of this experession (two for the +8 and two 
for the ) represent the four possible intersections of the 
panel leading or trailing edges and their extensions with the Mach 
line. A total of eight roots are obtained when both leading and 
trailing edges are considered. The scheme to identify the regions 
and limits on n takes these eight roots and arranges them so that 
they are in order according to their n value (from left to right). 
With each root is an identifier which states whether the root 
belongs to the leading or trailing edge. Roots to the left cf the 
left panel edge (y L ) are excluded as well as roots to the right of 
the right panel edge (y R ) . 

The integration limit determination procedure starts by 
establishing whether the Mach line crosses the left panel edge 
ahead, behind, or on the panel. If ahead of the panel, the first 
region, bounded by and the first intersection (which has to 
be with the leading edge) will be a Region 0 (no integration - no 
influence) . If on the panel, the first region will be a Region 
II. If behind the panel, the first region will be a Region I. If 
there are no intersections (roots) between y L and y R then the 
entire panel is of the same region type. The key region is Region 
II. From this region, the next region will be either 0 cr I 
depending whether the intersection is with the leading or trailing 
edge respectively. In going from one region tc another. Region II 
will always be involved. 

Figure 17 shows the logic block used to establish the limits. 
Once the procedure is started, at the left panel edge, it proceeds 
automatically with Yj, being the first limit and y R being the 
last. As an example of how the logic block results are used. 
Figure 17 also shows the regions and n limits of integration for 
a right going Mach line whose origin is positioned behind the left 
panel edge (R j denotes a root cf Equation (53)). 

One quick test which can be made to see if a panel has any 
influence on a control point before using the above logic is to 
check the value versus the most forward C value of the 
panel edge. If x. is upstream of this C value then the panel 
cannot influence the 1 control point regardless of the supersonic 
Mach number. This computation can save evaluating the sorting 
roots in many instances. 

Because the supersonic flow problem requires this added 
computation, which involves only the panel-control point 
relationship, and since each panel is covered by several Yj ,it is 
advantageous to integrate all the possible components'* of each 
spline function which share a particular panel-control point 
relationship. That is, the contributions to Aij from each of the 
twelve components of Yj in (Equation (25)) are not summed 
consecutively, but rather are added to as the program loops from 





REGION 0 REGION II REGION I 








one panel to another. This means that the total value of Aj . may 
not be available until all integrations over all the panels' 1 have 
been performed. Such an arrangement avoids setting the region 
limits more than once for any control point-panel relationship. 
This procedure offers a computational speed advantage over another 
possible arrangement which forms one complete A^j before going 
to the next (i.e., a complete Yj ). With this latter 
arrangement, the region limits must be computed several times for 
the same control-point panel relationship, due to the spline 
arrangement where functions occupy common panels (Figures 1C and 
12 ). 

This concludes the basic formulation of the numerical 
integration scheme used to evaluate A^* (termed the influence 
coefficient elements) . The next sectidn describes how these 
elements are used to match the boundary conditions for various 
flow problems. 
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5.4 Satisfaction of Boundary Conditions on Singularity Surfaces 

The preceeding sections showed the construction of the basic 
building block of the numerical scheme - tne vortex spline; This 
section shows how the strengths, Sj (multiplicative constants of 
each Y^), are determined by application of the boundary conditions 
on the singularity surface (region "s" of Equation (11))* Each 
element naturally satisfies the governing differential equation, 
since they were developed from the horseshoe vortex, and each 
element also satisfies all the continuity of vorticity required of 
the final solution (the purpose of the spline) , therefore the 
strengths of these elements are all independent. This 

independency means that every element can be separately adjusted 
such that it best satisfies the imposed boundary conditions. 
Since the splines have a limited number of degrees of freedom, 
they, in general, cannot exactly satisfy these boundary 

conditions. A compromise match must be accepted. The compromise 
used for the spline development was a least square error match. 
A solution based on this technique will yield the best overall fit 
to the boundary conditions. 

The square error involved in the approximation to the 
boundary conditions can be mathematically expressed as: 

e2= f f [\ c ('.")* w 8 k, ">] 2 dCdn <S4> 

8 

where is the square error, Wbc ( is the applied downwash 

boundary condition and Wg ( ? ( n ) is the downwash from the splines. 
Since it is not possible to perform the indicated surface 
integration for general cases, an approximation has to be made. 
This is dene by assuming the integral can be expressed as a 
weighted sum: 

* j 

M 2 

E 2 * E [yVV‘ W b t ,n )] ,55) 

1 ■ 1 

Here ( C , n j is one of M points on the surface (termed the 
control points/ and a, is an area associated with that pcint. 
(The choice of control point locations and their associated areas 
is discussed later.) Equation (16) can be used to express 

V V in Equation 9 ivin 9: 

M N _2 

E2 = 1 -1 [WfiC l " j?i AljSj 
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( 56 ) 


The rotation has been abreviated such that the subscript i 
denotes the value at ( ) . 

The square error in Equation (56) can be minimized with 
respect to each strength, Sj , since these strengths are 
independent. This minimization is accomplished by differentiating 


Equation (56) with respect to 
expression to zero: 


St 


and setting the resulting 


,1,2 M _ N 

si.. ' 2 E Kc, -EVhlV 


1=1 


1 


Ik 


(57) 


The variable, k, ranges from 1 tc N , therefore Equation (57) 
represents N algebraic equations which can be solved for the N 
strengths. A more useable form for Equation (57) is: 


M N 


M 


UTl A iw a i A ij s j = 52 A ik a i 

i- 1 j=i 


w 


1=1 


BC 


(58) 


1 


The entire set of Equations (58) can be written in matrix 
notation as: 

- - Mxl 


NxM MxM MxN 

[a] T N [a] 


NxM rp MxM f . 

- [a] M ( w bc) 


(59) 


ere f a! is the transpose of matrix [aJ with elements 


, anc£ J NiJ is a diagonal matrix with elements a . The 
, and ate vectors representing the unknown 1 strengths, 

, and the known boundary conditions, Wjjq . Solving Equation 
(59) for {S} gives: 

.-1 


l 


[ S | - rtAfNWl [a] T N[ W BC 

1 * H M 


(60a) 


Equation (6Qa) is valid whenever the number of control points 
is greater than or equal to the number of unknowns (M>N). If the 
number of control points is equal to the number of unknowns (M=N) 
then [A] is a square matrix and Equation (60a) reduces to the 
determinate form of: 


{s]-[ A ] W BC ] 


(60b) 


An abreviated notation is introduced here: 
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(61a) 


* T -= [a] T M 

A,[ A ] 



(61b) 


(61c) 


<61d) 


then 


, ^ T *1 «T . 
S * (A A) A 


(62 ) 


The symbol A T is the transpose of A weighted by M- 

When the special Mach line paneling is employed, the 
hyperbolic nature of the flow can be used to advantage by treating 
each region as a separate wing and solving the problem in steps. 
Figure 18 shows the arrangement of the panels with the associated 
[A] matrix partitioned into elements which represent the 
influence of subdivisions on themselves, A ; and elements which 
represent the influence of subdivisions cn other subdivisions, B. 
It is obvious from this figure that it is very inefficient to 
solve this problem with this one large matrix. Instead, the 
problem can be divided into three separate, smaller problems, 
corresponding to the three subdivisions (I, II, III). If this is 
done, the three equations to be solved in a manner analogous to 
Equation (59) are: 


AjSi 
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FIGURE 18 - PARTITIONING OF AIC MATRIX FOR MACH LINE PANELING 




(63b) 


A ra s ra * w m - w rn 


(63c) 


where 


w n - B n , i s i 


w m ‘ B ni,i s t **m,n s n 


(63d) 


(63e) 


The subscripts denote the subdivisions and the quantities Bj , 
and Wm are the known downwashes in each subdivision. The 
right hand sides of Equation (63b) and (63c) can be thought of as 
"effective" downwashes for the independently treated downstream 
subdivisions. 

Another computational consideration is the cost o ; f performing 
the operations indicated in Equation (60). A more direct method 
can be employed if Equation (56) is rewritten as: 


IN , 

£[ W Bc.V"®t “ £ a ij s j] 

1=1 1 J «1 


The above expression can be minimized directly using 
Householder's technique (Ref. 9). This procedure provides a great 
savings in computation time over the method represented by 
Equation (60) . 

The key feature of this least square error approach is to 
select the points ( C|, 0 |) and the areas, , such that Equation 
(55) does approximate Equation (54). with this goal, a study of 
various arrangements was undertaken. Fundamentally the study 
entailed subdividing the panels, which were used to form the 
splines. These subdivisions were formed by lines of constant 
percentage span and chord. Control points were located on these 
subdivisions using the area of the subdivision as the weighting 
factor, a| » in Equation (55). These points were located at the 
spanwise centroids of the subdivisions, midway between the leading 
and trailing edges. In the study the number of spanwise sections 
could be specified independently from the chordwise sections. 
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Figure 19 shows the minimum arrangements resulting frotr the 
study for supersonic and subsonic flow. By minimum, it is meant 
that the arrangement’s approximation to the integral (Equation 
(54)) was adequate over a wide variation of paneling. (Results 
presented later will show the accuracy of the arrangements fcr 
cases with as few as twelve panels and twenty-four control 
points.) The supersonic case required more control points because 
of the limited extent of influence in supersonic flow; 
consequently more points are needed to "sense” this influence. 

The study also revealed that difficulties are encountered in 
the vicinity of subsonic leading edges and at all tips due to the 
singular nature of the vorticity distribution in these regions. 
Since the splines do not have the proper asymptotic form, the 
downwash has an improper singular behavior. This downwash 
departure tends to have a strong influence on the error 
minimization procedure if ccntrol points come too close to the 
edges involved. So, another study was made to determine suitable 
locations for control points in these regions. A review of other 
lifting surface schemes showed that their respective two- 
dimensional formulations could give the proper lift on a flat 
plate at angle of attack in subsonic flow using only one element 
(e.g., N- 1 in Equation (16)). in fact, this is a criterion that 
can be used to establish the proper control point location fcr a 
determinate numerical formulation (where the number of unknown 
strengths is equal to the number of control points). This two 
dimensional criterion can be shown to be valid for the vortex 
lattice (Fef. 10) and the constant pressure panel schemes. Using 
this criterion, the proper control point location tor the spline 
was computed. The one element spline representation for the lift 
on a flat plate is shown in Figure 20. The integrated load for 
the spline was set equal to the linear theory integrated load. 
The location was then found where the downwash from the spline is 
equal to the downwash on the flat plate. This location is at the 
center of the panel. (This result can be confirmed using Eguation 
(18) with the distribution T(C) * ( x 2 " C )/( x 2 “ x i) )• Starting 
with this result, it was discovered, by numerical investigation, 
that points on or downstream of the center chord line of the 
leading edge panel consistently produced accurate results. Since 
the components of the spanwise functions begin and end at the mid 
panel span, (Figure 10), the criterion in this direction is that 
the control points near an edge be on or inboard of the three 
quarter panel span (corresponding to the center of the half 
panel) . 

Using the above criteria, the general paneling arrangements 
shown in Figure 8 could be applied directly without any special 
considerations for subsonic flow, since the control point 
configuration satisfies the conditions outlined. For supersonic 
flow, there is no problem spanwise; chordwise, with a supersonic 
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FIGURE 19 - CONTROL POINT CONFIGURATIONS ON A SINGLE PANEL 







Leading edge, there is also no problea. But, if the leading edge 
is subsonic, the two forward control points oust be excluded from 
the ainimization process. With this siople exclusion rule, 
consistently accurate results were obtained over a wide range of 
test cases. 

The benefit of this least sguare error scheme can be aore 
fully realized for cambered surfaces where the nuaber of functions 
needed to model a problea can be small compared to other schemes. 
This is because relatively few functions can model rather rapid 
camber variations and the number of functions is not tied one to 
one to the nuaber of control points. With the least square error 
formulation, any number of control points per panel can be 
introduced to accurately incorporate the camber variation into the 
problem and at the same time, the number of functions can be held 
to a minimum. As pointed out earlier, the configurations shown in 
Figure 19 represented the minimum arrangement for consistent 
accuracy, any larger nuaber of control points could also be used 
to obtain better camber definition, (subject to the edge 
constraints outlined above) . 
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5 . 5 Computational Procedure 


Figure 21 shows a flow diagram of a method to compute results 
using the vortex spline scheme. Once the strengths, S, have been 
computed, the remaining parameters can be determined (e.g., lift, 
moment, pressure distribution, etc.) . Since the splines provide 
a continuous distribution over the surface, the pressure and bound 
vorticity can be evaluated directly at any point by determining 
the aggregate influence of the contributing functions. This 
continuity also makes it possible to compute the shed vorticity, 
a , by integrating and differentiating the basic distribution for 
y as indicated by Equation (23) . 

The lift and load distribution can be computed either by 
actually integrating the basic functions or by performing a sum on 
the pressures, similar to what was done for the least square 
downwash error integral formulation using the same control pcints 
and areas. The moment can likewise be computed by integrating or 
summing; however, the accuracy of this calculation using the 
summing technique is less than that for the lift calculation. 
Nevertheless, the summing technique is sufficiently accurate when 
several panels are employed and was the method used to compute 
both the lifts and moments presented in the next section. 

In addition to the standard aerodynamic parameters, the mean 
square error can also be computed. This mean is the parameter, E? 
in Equation (55) divided by the total area (sum of the ). The 
mean square error (with appropriate control point exclusion) is an 
excellent indicator of the overall accuracy of the solution. The 
square root of this value, the root mean square downwash error, 
was found to be the same order of magnitude as the root mean 
square pressure error for the cases examined. This BUS error 
parameter supplies the user with an a posterior judgment of the 
validity of a flow model without the need to consult other data. 
It also opens the possibility of formulating a scheme which would 
automatically adjust the paneling and rerun a case if sufficient 
accuracy were not obtained from the initial run. 

The next section presents the test cases that were used to 
evaluate the scheme. 
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Mach Line Paneling Loop 


Solution Continues 
for Additional Regions 


“Effective" 
Downwashes Used 
For Downstream 
Regions In Super- 
sonic Flow 
(See egn.'i 63b 
and c) 


Geometry 


o Areas 

o Control Point Locations 
o Panel Descriptions 


. AIC Matrix ? : 

Generator 

o A A B Downwash at Each Control Point Due to Each Function 

(A - Aerodynamic Influence Coefficient Matrix)* 


Least Square 
Error AIC Matrix 
Formula tor 


o a” T A 


Least Square 
Error Downwash 
Matrix 


Solution- 
Strengths of 
Singularities 


o A^fW) 


-1 


O {8J-[a T a] A T (W) (eqn. 62) 


0 (W) - B{8) Downwash at Points 

eqn.'s 63 d In Downstream Regions 
and a) ~T~ 


I Bound Vortlclty /vV 
Calculations '■ ’ 


o Total C. , C M & Wing Bending Moment} Momt. A Load Dlst 
EXIT (After Additional Points) 




Yes 


Additional 

Regions? 


JESL 


Error Analysis 

} oErrors Between Solution and True Downwashes 
Additional points 
at which the vortlclty Is desired 

Exit 


* The Influence of the functions on points outside the region of Investigation 
(Matrix B) are also computed. These points are used In the supersonic case to 
compute the influence of one region (upstream) on another (downstream). This 
technique allows the solution to be obtained by marching downstream, (see flg.18.) 

o- Computed Quantity 

(**j Denotes Weighted Matrix 


FIGURE 21 - CONCEPTUAL FLOWCHART (FIRST ORDER) 












5.6 Numerical Besults 


Various configurations were analyzed to determine paneling 
requirements and convergence to a known answer. The spline 
function results were compared with other lifting surface schemes 
and exact, linearized flow solutions. Figures 22 through 28 show 
the results. 

Figure 22 is a comparison of the spline scheme with another 
lifting surface solution (Ref. 10) for a circular wing in 
incompressible flow. The comparison shows good agreement not only 
for the coefficients but also for the detailed loading. Note how 
well the quadratic edged panels can match the exact planform. 

Figures 23 and 24 show comparisons of the spline scheme with 
a pressure mode scheme (Ref. 2) for compressible, subsonic flow. 
Again agreement is very good. Note the rounding of the planform 
representation at the root for both cases. This eliminated 
downwash singularities without introducing any special singular 
forms for y . 

Figure 25 ' shows a comparison of the spline scheme, using 
special Hach line paneling, with the exact, conical flow result 
(Bef. 11). The planform paneling shows a rounding at the root. 
This rounding produced a slight departure of the pressure from the 
exact result in the vicinity of the rounding (not noticeable in 
the chordwise plot of the pressure) . 

Figure 26 'shows two comparisons of the spline scheme using 
geometric paneling for supersonic flow: one with the exact 
solution and one with a constant pressure panel scheme using the 
same paneling. Even for this paneling, which ignores the special 
Mach lines, the spline scheme gives accurate coefficient results 
with a reasonable comparison tc the pressure. 

Figure 27 shows another comparison of the spline scheme with 
exact conical flow results (Ref. 11). This planform has a 
subsonic leading edge. Again note the rounding of the planform. 
This rounding causes a local departure of the pressure from the 
exact result as shown. The overall agreement is good, even for 
the relatively few panels used. 

» *‘ > 

Figure 28 shows the numerical stability of the formulation 
for various paneling schemes. In all cases the pressure 
distributions appear to be as accurate as could be expected for a 


65 





k a= 0.0995 RADIANS PRESSURE MODE (REF. 2) 

;.'M ; =0.12 VORTEX SPLINE - 



(N„ is the neutral point expressed 
in terms of fraction chord aft of 
root chord leading edge.) 


Root Mean Square Downwash Error - 2.15 X 10~^ 


■FIGURE 23 - UNTAFERED SWEPT 'WING 
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a = 0.199 RADIANS PRESSURE MODE IREF. 2) 

M = 0.15 VORTEX SPLINE 



(N 0 is the neutral point and is • ■. i • 
expressed in terms of fraction mean 
chord aft of root chord leading edge.) 


Root Mean Square Downwash Error = 1.81 X 10 - ^ 


FIGURE 24 


TAPERED WING WITH ZERO SWEEP OF QUARTER CHORD 




Root Mean Square Downwash Error =4.63 X 10 


FIGURE 25 - SQUARE WING (SPECIAL MACH LINE PANELING) 




Ci * 2.07 

CONSTANT PRESSURE PANELS „.357 


Root Mean Square Down wash Error - 1.79 X 10' 2 


FIGURE 26 - SQUARE WING (GEOMETRIC PANELING) 
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o=1 RADIAN 

M = VI 


~EXACT(REf.11,PG.157l 
— VORTEX SPLINE 



fraction of mom chord ift of opox.) 

_2 

Root Mean Square Downwash Errors 5*00 X 10 • 


FIGURE 27 - DELTA WING 





Root Mean Square Downwash Error =1 . 02 x 10 



Root Mean Square Downwash Error = 1.02x10 


-1 


FIGURE 28 - VARIOUS PANELINGS OF A CIRCULAR WING 



given paneling using a linear chordwise representation. The 
stability of the solution does not depend on paneling, and the 
solution converges to the exact results as the number of panels is 
increased. 
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6.0 SHELL REPRESENTATION USING, THE SPLINE, SCHEME 


6.1 General Approach 


This section presents a scheme which represents the 
interference effects on wing-body combinations. The basis for 
this representation is presented in Beference 1 where it is 
referred to as an interference shell. The shell amounts to a 
cylindrical surface with singularities distributed over it, much 
the same as for the wing. The interference shell, however, is not 
a impermeable surface with flow tangency boundary conditions. The 
purpose of the shell is to account for the mutual interaction 
between the body and the wing. The shell boundary condition is 
that the shell singularities induce a velocity which cancels the* 
velocities on the shell that are induced by the wing. The body 
contained inside the shell, will however, induce a flow through 
the shell which is not canceled. The enclosed body then appears 
to be in an undisturbed free stream flow, therefore its analysis 
can be treated separately using some appropriate technique (e.g.. 
Ref. 12) . 

The aim of this shell study was to provide an improvement 
over the constant pressure panel shell (fief. 1). In so much as 
the vortex spline arrangement provided much improvement over the 
constant pressure panel scheme for the wing, it' was felt a similar 
development for the interference shell would be appropriate. 
There was, however, an unresolved question concerning the 
interference shell - namely whether the elementary horseshoe 
vortex was the appropriate singularity needed to model the body- 
wing interaction effect. The question arose because bodies in 
potential flow do not have a wake (i.e., shed vorticity) and 
because the elementary horseshoe vortex does have a wake. This 
difference between the flow cn a wing and the flow on a body can 
be attributed to the roundness of the aft end of a body as 
compared to the sharpness of the trailing edge of a wing 
(suggesting no Kutta condition on a body). If the Kutta condition 
is not used, a vortex solution is not unique (fief. 5, pg . 57). 
This means that a shell representation using the vortex spline 
which needs the Kutta condition might not properly model the flow 
on a body. Decause of this possibility, a study was initiated 
whereby several elementary singularities were developed which 
could be utilized in a spline technique. These singularities were 
the elementary horseshoe vortex, the source, and the doublet. 
Unfortunately, none of the analyses was developed to the point of 
yielding numerical results. The next section shows the analysis 
which was performed for the various singularities as well as a 
detailed description of a scheme using the elementary horseshoe 
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vortex in an influence coefficient approach. The scheme described 
is compatible with the scheme already presented for wings. 



6.2 Interference Shell Singularities 


Figure 29 shows the cylindrical coordinate system used to 
analyze the interference shell. Since the doublet and elementary 
horseshoe vortex can be derived from the source (fief. 5, pp. 28 
and 87) , the analysis will begin with the source singularity for 
subsonic flow. 


The elementary source potential , $ , 

given as (fief. 5, pg. 87): 8 

*«,(* , r ,e ) = 1 ? 


for subsonic flow is 


4 ir 


m 




(65) 


where p^ = l-M^ • in this expression, the source has been placed 
at the origin of the coordinate system. If however, the source is 
located at some other point ( x , r ,6 ) , then the potential 
becomes e dependent: ° 0 


*s (*> r . e , x o > r o » 0 o ) = " "" — ~ (66 ) 

4 * ^(x-x 0 ) 2 + B 2 (r 2 +r2_2r r Q co8 ( 6 - e o )) 

The above result can be better understood from Figure 30. 

The kernel function corresponding to that used for the planar 
vortex panels can be derived by differentiating (6fc) with respect 
to the local surface normal coordinate (i.e. , 9$ g /3r |r»r Q is the 
downwash on a cylinder of revolution). In order to obtain the 
influence coefficient, Ajj this velocity kernel must be 

multiplied by the proper strength variation chordwise and 
circumferentially (the latter variation corresponding to the 

spanwise variation for planar panels). After multiplication, the 
function must be integrated over the domain of the singularity. 

The interference shell can be paneled in a way entirely 
analogous to the geometric paneling scheme used for the wing 
(Figure 31). First the cylinder can be divided streamwise by 

parallel rays on the surface. Then the cylinder can be further 

divided by parallel planes normal to the axis of the cylinder. 
These cuts define a group of panels which are rectangular in the 
x-0 coordinates. Spline functions for the various 
singularities can be defined using these panels. 

Since the different singularities require different degrees 
of continuity, the analysis from this point or. will be limited tc 
the elementary horseshoe vortex using the continuity conditions 
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FIGURE 30 - AXIS SYSTEM AND DEFINITION OF VARIABLES 




FIGURE 31 - PANELING OF AN INTERFERENCE SHELL 


e 



FIGURE 32 - THREE PANEL SPLINE FUNCTIONS ON INTERFERENCE SHELL 
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already defined for planar panels. The analysis will be further 
restricted to circular cylinders having equal arc length panels. 


The first step in the shell analysis is the derivation of the 
elementary horseshoe vortex on a circular cylinder. This 
singularity is obtained by integrating the doublet potential 
function (Ref. 5, pg. 87) from some initial point Xo to » 
(downstream) . This integration corresponds to a line of doublets 
from the point Xo to « . The potential ot the doublet with its 
axis in the ? Q direction (Figure 30) is derived by 

differentiating the source potential (66) with respect to r c ' 




[(x-x„) 2 4S 2 (r 2 + r| -2rr 0 cos( 8 . 9<> )] /2 (6?) 


The potential of the elementary horseshoe vortex is then: 

00 

♦ v (*» r .•» *o> r o » # o ) - j ♦(*,*,•, * •© )**; 

*o 

r 0 -rco»(>-«„> f. »-»« 1 (68 

4. (r 2 4r 2 -Irr 0 co. ( ,-. 0 l|( /,x-, o ) 2 , , 2 ,r 2 * r 2 . 2rr co.(. -.)> 


For this analysis, a spline arrangement spanning three 
circumferential panels was used (Figure 32) . The circumferential 
variation can be derived using three quadratics (one on each 
panel) and a set of constraint requirements. On the first panel 
(Fig ure 33) , the requirements are that the function must vanish in 
value and slope at the left edge. On- the second panel the 
requirements are that the function be continuous in value and 
slope at the left edge. On the third panel the requirements are 
that the function be continuous in value and slope at the left 
edge, and vanish in value and slope at the right edge. These 
eight requirements, coupled with the additional specif ications of 
unit peak amplitude, completely determine the nine coefficients of 
the three quadratics composing a single spline. The spline 
distribution of bound vorticity (strength of elementary horseshoe 
vortex) can then be written as: 

, . « ^ I m 

Y j x,fl 'liL )0 ^ (69) 

where x! is the same as that defined by Equations (26) and (27) , 
and the 9 QT for equal arc length panels are: 
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_ 1 2-2 

Vt 9 , 


( 70 ) 


Q ) 2 ' 


— 4 — e e 2 

3 3 2 3 2 


(71) 


Q 3 ^ _2 4 _~ + 2-2 

°J 3 3 6 3 3 ®3 


(72) 


The 9 s are local coordinates for each panel ranging 
0 at the left edge to 1 at the right edge. 


from 


The influence coefficient kernel is found by differentiating 
Equation (68) with respect to r (to find the normal velocity 
kernel) and evaluating this resulting expression on the surface of 
the cylinder ( r equal r n ) • 


K( V e 1 »*o»*o» T o 1 




i ~ 2 co she 
2t£ (1-cosAfl ) 


f x l“ x o-l 

l 1+_ R J" 


( x j~ * 0 )(i-CosA6) 


2R3 


] 


(73) 


where Ae is« r 0 O , R is ^ x l“ x o^ + 2 flV (i -cosAO), 
"i” denotes that the field point ( r , 0 , r ) is 
surface. 


and the subscript 
on the cylinder 


The influence coefficient is then: 




(x, e ) KUpe^x ,0)dxrd0 


(74) 


The integration can be carried out in exactly the same manner 
as for the planar panels performing the x integration first then 
using the numerical scheme, previously described, to integrate in 
the 0 direction. This same numerical integration scheme can be 
used since, the cylindrical kernel has the same singular properties 
as the planar kernel. 
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For .supersonic flow, the differences in the cylindrical 
analysis are the same as outlined for the planar analysis 
basically the kernel and limits of integration change. 


The supersonic source is (E?ef. 5, pg. 87) : 

9 2« / 2 « — o Z~~~ 

V(x-* 0 ) - ^(r 2 4r o 2 -2rr 0 co»(»-«^) 

The doublet in the ? 0 direction is defined by ♦<! ■ : 


(75) 


3r_ 


A Iw r m w - Mi tf 2 (r 0 .rco8(«-e ft y) 

♦*<*»*»• .*o» r o» # o» — r :: — r*— 


2w [(x-Jf 0 ) 2 -i 2 (r 2 4r 2 - 2rr 0 cos(#-e 0 )] 

and the horseshoe vortex is defined by 


H 


(76) 


-J ** 


- ^t 2 (r 2 + r 2 - 2rrco«(» 


dx. 


(77) 


r t i JL r o- rco<, < g - > ») 

% »*o» r o*'ol ». « • 


x-x. 


2" (r*4r2. 2ProCO ,„. v , y, , 2. ,2 (r 2 <r2 . 


2 r r coi ( • -d^> 


In the above 0 2 is defined as M 2 - 1 ( 8 2 is always positive 

whether supersonic or subsonic) . 


The kernel needed for the. influence coefficient is derived 
from (77) by differentiating with respect to r and evaluating 
the result on the surface of the cylinder (r equal r ) : 

O v 


^ » 8j » 


X O. 0 n> r c 




2 (1-cos A0) 


(Xj-x 0 ) 

R 


- 8 


2 (xj-x^Ki-eoBAe) 


2R 3 


(78) 


where ao is e t - e G and R is , -x© f~ 2 0 2 r c 2 f 1 - cot A6) 

The limits of integration for the supersonic integral are 
found by determining the intersection of the forward Mach cone, 
emanating from the control point ( Xj , 8^ ) , with the cylinder. 
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The equation of the Mach cone evaluated on the surface of the 
cylinder is: 


For the simple panels used 
are found by replacing x with x (x of the leading edge) and 
* T £ (x of the trailing edge). Equation (79) then gives the 8 
values of the intersections which can be used in the integration 
limit scheme previously described to establish the region (0, I, 
or II) and set the limits for the 6 integration. 

Once the influence coefficients are computed, the boundary 
conditions are satisfied and remaining analysis is performed using 
the same logic as was used for the planar development. The 

boundary conditions of course are different for the interference 
sell as mentioned before and are more fully covered in Beference 
1. However, this difference does not affect the mechanics of 
obtaining a solution, it only changes the layout of the overall 
influence coefficient matrix in that the body contained within the 
shell is assumed to have no influence on the shell. 

The above analysis could be repeated using doublets, but 
since the equivalent doublet strength is the integral of the bound 
vortex strength in the x direction (fief- 5, pg. 135), such a 

formulation would require quadratic distributions in both X and Y. 
This arrangement could be made using the relations given by 

Equations (70) through (72) for both the X and 8 variations. An 
influence coefficient scheme could then be formulated which would 
not shed vorticity. This scheme could then be compared with the 
vortex scheme to determine the merits of each for body 

representation. The scope of the present development did not 
include such a comparison. 


2a 2 r* [l-cos(8|.8)] 


(79) 


in the analysis, the intersections 
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7.0 SECOND ORDER SOLUTIONS USING THE VORTEX SPLINE 


In the introduction it was pointed out that one of the 
purposes for developing a new aerodynamic technique was to have a 
scheme which would calculate those lateral-directional stability 
derivatives which cannot be handled by th 2 current formulations 
which treat arbitrary wing-body combinations. The previous 
sections have shown the development of a scheme which was felt 
would provide the necessary improvements to meet the lateral- 
directional requirements. The course of development was to be an 
initial treatment of planar wings (Ref. 13), then a treatment of 
bodies in the presence of wings (the interference shell) , and 
finally the treatment of these wings and bodies and their 
combinations in interaction flow. The initial interaction flow 
studied was for angle of attack-angle of sideslip. 


The basis for the second-order flow analysis is given in 
Reference 4. In essence, the development assumes that the exact 
nonlinear, potential flow equation can be expanded in terms of 
small geometric parameters which define the problem (e.g., angle 
of attack, angle of sideslip, wing thickness, etc.) . When this 
expansion is carried out, an ordered set of linear equations is 
obtained. Each of these equations has its own boundary conditions 
which are determined by a similar expansion of the exact boundary 
conditions. The higher order equations (those involving powers 
and products of the small parameters) have non- homogeneous 
differential equations while the first order equations have the 
familiar linearized, homogeneous differential equation, e.g.. 


(l-M 2 ) 





f(x, y* z) 


(8 C) 


<>-«’> * &-• ,10 ’ 

where the subscripts "I” and " 2" denote the first and second order 
terms from the potential flow expansion. The expansion assures 
that F(x,y,z) is only a function of ^ and its derivatives. 

The forms of these results (Equations (80) and (10)) are 
standard ones arising from expansions so that the technique used 
to evaluate results for tnese equations is general. This point is 
being made because the expansion in Reference 4 did not properly 
account for all the terms present and therefore is not valid. 
Although the result is not correct, the forms of the equations are 
indicative of the types of numerics needed. The work done in 
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implementing the technique is therefore valid and would be 
applicable when the proper expansion is found. Because of this, 
the technique used to solve the class of problems represented by 
Equation (80) is presented, but the presentation will be somewhat 
general in nature. Detailed extension of the scheme to second 
order wing-body configurations was of course not done. 

The solution to Equation (80) consists of two parts - a 
homogeneous part, , and a particular part, <£ 2 p . Since 

Equation 10 can be added to Equation (80) in any proportion 
without changing the problem, an assumption is made that a v 2 p can 
be constructed which will satisfy the wake boundary conditions. 
This leaves homogeneous boundary conditions in the wake for the 
homogeneous solution $2H • This means that the same homogeneous 
formulation, (i.e., vortex spline) can be used for both first and 
second order homogeneous solutions. The particular solution, 
$ 2 p , provides a basic velocity kernel in exactly the same way 
the horseshoe vortex (homogeneous solution) provided the velocity 
kernel for the wing and the interference shell (Equation (68)). 
Following the procedure outlined for the shell, this potential 
function for ^ 2 p is differentiated with respect to the surface 
normal at the control point. The resulting differentiated 
quantity is the normal velocity kernel, analogous to those 
associated with the first order problems for wings and 
interference shells (Equations (12) and (73), respectively). 

Unlike the first order problem the particular solution does 
not depend on matching the boundary conditions. Bather, the 
strengths of the distributions and the distributions themselves 
are known directly and no solving is involved. This is because 
F(x,y,z) is only a function of the first order problem, whose 
strength is determined initially in the solution to that problem. 
Specifically then, the velocity kernel determined from is 

used to compute an additional influence on the boundary conditions 
and not used as a means to satisfy them. The downwashes computed 
from the particular kernel and those obtained from the boundary 
condition expansion are combined to provide a new set of boundary 
conditions for the second order homogeneous solution. The 
homogeneous solution, <Pin , is then used to satisfy the combined 
boundary conditions. The net result is that the combination of 
^ 2 p and ^2H satisfies the boundary conditions for the second 
order problem which were determined by the initial expansion. 

Since the second order homogeneous differential equation is 
identical to the first order equation, the same basic solution 
technique can be used (i.e. the vortex spline). The second order 
homogeneous solution, of course, has different boundary conditions 
and correspondingly different strengths from that of the first 
order. Also, depending on which equation is involved, there can 
be different symmetry requirements for the solution (e.g. angle of 
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attack -angle of sideslip problem assuming zero sideslip is 
antisymmetric while the first order problem, angle of attack, is 
symm et ric) . 

Figure 34 shows how a seccnd order problem could be solved 
using the vortex spline scheme. Basically the central block of 
logic is the homogeneous boundary condition solver which is the 
vortex spline scheme (Figure 21). The second order particular 
solution section of the logic blcck provides the downwashes which 
must be combined with the specified boundary condition (from the 
expansion) . These downwashes then represent a new set of boundary 
conditions for the homogeneous solution ( ^ 2 h ) which is solved 
in the exact same manner as the first order problem. The 
quantities necessary to obtain the pressures are then computed 
(this is not always just the bound component of vorticity as 
pointed out in Reference 4) and the stability derivatives are 
obtained. 

The advantage of the vortex spline for solving these second 
order problems can be seen from a study of the characteristics of 
the forcing functions F(x # y,z) (Ref. 4). These functions 
generally involve derivatives of the first order solution. A 
solution to the first order problem must, therefore, be capable of 
providing sufficient derivative resolution that the second order 
problem can be solved accurately. The splines provide this 
capability as can be seen from Figures 22 through 28. 
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KHiJ^aointUm strengths Computed from First Order Solution 1.e.eq.(62) 

o Downwashes at Control Points Due to Particular 
I 'Solution 

CThls section is the saa» as for the first order 
solution except that the boundary conditions 
are established by the particular solution 
and not by the geometry and angle of 
attack of the model.) 



Error Analysis 



o Average Error t Average Square Error 

I Circulations from Homogeneous I Particular 
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*The Influence of the functions on points outside the region of 
Investigation (Matrix B) are also coeluted. These points are 
used In the supersonic case to compute the influence of one 
region (upstream) on another (downstream). This technique allows 
the solution to be obtained by marching downstream. 
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FIGURE 34 - CONCEPTUAL FLOWCHART (SECOND ORDER) 








8.0 CONCLUDING REMARKS 


The experience gained with the vortex spline scheme has 
indicated the following: 


1 . 

The 

solution 

is 

not strongly 

dependent 

on paneling 


a rra 

nge me n t ; 

i t 

works well 

wit h 

uniform or 

nonuniform 

panel 


spac 

ing . 







2. 

The 

scheme 

works 

equally 

we 1 1 

for both 

subsonic 

and 


supersonic flow. 

J. Special Mach line paneling produces better pressure 

resolution in supersonic flow than geometric paneling. 
However, in all cases the geometric paneling produced 

accurate values for lift and moment. , 

4. The absence of downwash singularities in the wake (except 

behind planform tips) implies that tails can be paneled 
independently from wings. 

5. Computation times are small, requiring only a few seconds 

(CDC 6600) for simple planforos and less than 2 minutes for 
finely paneled planforms. 

6. Accurate results can be obtained with sparse paneling. 

7. The extension of the scheme to the representation of 

arbitrary configuration surfaces appears entirely feasible. 

8. The numerical stability of the scheme provides flexibility in 
the application of the method. Such stability is a 
fundamental requirement for the computation of downwash 
(camber slope) given pressure (this type of problem solving 
is termed the design mode) . 


8. The scheme is self-checking in that the user can judge the 
accuracy of the numerical representation by simply examining 
the root mean square error in the downwash. (This rms value 
seems to be the same order of magnitude as the rms" error for 
the pressure distribution for the cases tested sc far.) 

Future plans for the scheme should include investigation of 
an automatic paneling technique wherein the user would only be 
required to input the vehicle shape and specify an acceptable rms 
error. The program would then automatically panel the 
configuration, solve the problem, examine the rms error, and 
repanel based on the error. This cycle would continue until the 
specified rms error was achieved. such a development would lead 
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to' a self-checking scheme that provides results of- a predetermined 
level of numerical accuracy. 
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10.0 APPENDIX 


Doublet Characteristic Box Method 


This method is an alternate approach to the vortex spline in 
supersonic flow. The method shoved promise but further 

development was halted when it became evident that a great deal of 
time would be needed to calculate all the integrals involved for 
the various building blocks, as well as developing a systematic 
approach to identify all the "special" Mach lines. There was, 
however, substantial work done on the fundamentals. This effort 
is presented below. 

The first step in the present work was to identify the 
"special" Mach lines emanating from planform edge breaks. These 
lines are used to divide the planform into several different 
regions. Each of these regions is divided into a network of 
panels bounded by Mach lines as shown in Figure 35. 

The basic building blocks comprising the method were 
identified. These blocks, along with their controlling 
parameters, are shown in Figure 36. The building blocks include 
(see Figures 35 and 36) : 

I Regular - Bounded by Mach Lines 

II Two Subsonic Leading Edges 

III One Subsonic Leading Edge 

IV Supersonic Leading Edge 

V Subsonic Trailing Edge 

VI Wake 

A planform is represented by a superposition of these 
building blocks. The singularity distribution on each panel is 
assumed to be a function of the Mach line coordinates (r, s) 
Figure 37. The type of singularity distribution that was explored 
initially was a doublet distribution of the form 

= a + br + cr 2 + ds + ea 2 + fra + gr 2 s + hra 2 + kr 2 a 2 ( 81 ) 


where &+ is the velocity potential jump across a panel and 
a,b, ... , k are unknown coefficients. This distribution provides 
vorticity components that vary linearly and quadrat ically with r 
and s as shown in Figures 38 and 39, and insure smooth derivatives 
of the resulting first order solutions. 
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FIGURE 36 - BUILDING BLOCKS 
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FIGURE 37 - CHARACTERISTIC COORDINATES 


"Special" Mach Lines 



FIGURE 38 - INTERIOR REGION PANELING 



vorticity component parallel to r 
vorticity component parallel to s 

varies linearly with r 
Y s varies quadratically with s 

Y r varies linearly with s 

Y varies quadratically with r 


FIGURE 39 - BASIC VORTICITY DISTRIBUTION 
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The unknown coefficients are found in the following ways: the 
solution starts with a leading edge panel (II or IV Figure 36) for 
which hi is known (0 at leading edge); the solution then proceeds 
downstream by matching at common edges. These matching conditions 
consist of requiring hi , Y r * and Y s (Figure 39) to be 
continuous across the panel edges. This leaves one free parameter 
remaining for each panel which is determined by the downwash 
boundary conditions. The basic method of analysis is shown in 
Figure 40 for the "Regular" panel. The procedure consists of 
matching A<t> and the two derivatives, 3A<J>/dr and dhi/ja* along the 
two upstream edges. This leaves only the parameter k as unknown, 
k is then determined by satisfying the downwash boundary condition 
at a control point placed on the panel. 

Along the network boundaries (either "special" Mach lines or 
planform edges), some of the matching conditions (i.e., A<f>, BA<f> /3r 
and aAtf/38 ) are relaxed. They are replaced by additional downwash 
boundary conditions to provide the freedom needed to account for 
loading discontinuities, Kutta conditions or other constraints. 
Although the above procedure provides a determinate set of 
equations, a minimum error technique might also be incorporated to 
enable multiple boundary points cn the panels. 

The influence coefficients relate potential jump, a$ , to the 
downwash , a , through the planar relationship (the basic doublet 
representation - Ref. 14). 


(82) 


Here ft is the Mach number and [!„ is the free stream 
velocity. 

Thus, the unknown coefficients are determined by combinations 

of: 

• Downwash (camber) 

• Matching - hi , , h Cn 

3r as F 

• Minimum Error 

An investigation of the various building blocks and the 
matching conditions provides the results shown in Figure 36. This 
figure illustrates the matching conditions applied and the number 
of resulting free parameters, any of which might be relaxed. Note 
that a wake panel is completely specified (no free parameters) by 
the upstream conditions. 


2vu 00 a 


* 2 __ [( hi (r 1 ,s 1 )dr 1 da 
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A# and 


continuous 



3A» 
3 r 


A# = a + br + cr + ds + es^ + frs + 


gr 


Y s * 4^ ■ b + 2 cr + fs + 2 grs + hs 2 + 2 krs 2 
<7 r 

Y r - = d + 2 es + fr + gr 2 + 2 hrs + 2 kr 2 s 

o P 


1) At r = 0, given, 

2 

a + ds + es matches the upstream quadratic function 
a,d,e determined 

2) At s = 0, M given b and c determined as in 1) 

3) At r = 0, given b + fs + hs 2 f,h determined as in 1 ) 

o X* 

4) At s = 0, Mi. given, d + fr + gr 2 g determined as in 1) 

38 

5) Jc unknowns determined from the downwash constraint 

a~//A4K ( r , r-j ,s,s^ ) d^ ds ] 

where K (r,r^,s,s-|) is the kernel from equation 82. 


FIGURE 40 - "REGULAR" PANEL MATCHING 
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With the sideslip problem in view, the following basic steps 
were used to make up the doublet characteristic box method. 

1) The functional form for the potential jump was selected 
(i.e.. Equation (81)). 


2) The boundary conditions are used to determine the parameters 
of the doublet distribution through combinations of known 
downwash (camber) , matching conditions and methods of minimum 
error. 

/ 

3) a) Pressure (and ensuing load) can be found directly as the 

spatial derivative (in the free stream direction) of the 
potential jump. 


b) Second order sideslip calculations are made possible by 
the representations of the velocity potential through 
the improved spatial derivatives of the first order 
potential function. 

The complexity of the doublet characteristic box method 
hinges on the downwash- potential jump relationship which must be 
integrated analytically, numerically, or a combination of the two. 
Accordingly, the basic downwash- potential jump relationship might 
be cast into a more suitable form than presented by Equation (82) 
- such a form might carry the double differentiation under the 
integral sign, for example. 


An investigation of the various building blocks and the 
matching conditions provides the results shown in Fiaure 36. 
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